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>>: 

^ ' Design and reconstruction of 2d and 3d photonic structures are usually 

P-i' 

carried out by forward simulations combined with optimization or intuition. 
Reconstruction by means of layer-stripping has been applied in seismic pro- 
^ I cessing as well as in design and characterization of Id photonic structures 

^ : such as fiber Bragg gratings. Layer-stripping is based on causality, where the 

O ■ earliest scattered light is used to recover the structure layer-by-layer. 

^ ■ Our set-up is a 2d layered nonmagnetic structure probed by plane polarized 

harmonic waves entering normal to the layers. It is assumed that the dielectric 
permittivity in each layer only varies orthogonal to the polarization. Based on 
^, , obtained reflectance data covering a suitable frequency interval, time-localized 

'. pulse data are synthesized and applied to reconstruct the refractive index pro- 



file in the leftmost layer by identifying the local, time-domain Fresnel reflection 
at each point. Once the first layer is known, its impact on the reflectance data 
is stripped off, and the procedure repeated for the next layer. 

Through numerical simulations it will be demonstrated that it is possible 
to reconstruct structures consisting of several layers. The impact of evanescent 
modes and limited bandwidth is discussed, 
(c) 2013 Optical Society of America 

OCIS codes: 050.5298, 290.3200. 
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I. Introduction 

Photonic crystals liave been an exciting field of research since Eli Yablonovitch and Sajeev 
John published their papers in 1987 [HE]- Photonic crystals are periodic structures designed 
to affect the propagation of electromagnetic waves [3] . The usual way of constructing such 
optical components is to carry out numerical simulations of electromagnetic waves hitting 
and being scattered off a trial design. The design is then changed by optimization or intuition 
until the structure has the desired properties. In this paper we use a layer-stripping procedure 
to show that it is possible to reconstruct a photonic structure from a set of reflectance data 
based on harmonic waves. In other words, we will look at how the inverse problem can be 
solved. The ultimate goal is to be able to reconstruct the structures from the observations 
and use this as a supplementary tool for the design of structures with desired features. 

The idea of layer-stripping is to reconstruct the properties of a medium from scattered data 
originating from an emitted wave pulse on the boundary of the domain. Layer-stripping is 
based on causality, where the earliest scattered wave from each depth layer is used to recover 
the structure layer-by-layer. The method was first used in solving inverse acoustic scattering 
problems for seismic data; overviews can be found in [HIS]. Later, the method has also been 
applied to the design and characterization of one- dimensional photonic crystal structures 
such as fiber Bragg gratings [6] and multimode structures [7]. 

In the present simulation study, we are considering a two-dimensional structure, layered 
normal to the z-direction, with periodically varying refractive index in the x-direction while 
being constant in the y-direction. The structure is probed by plane polarized harmonic 
electromagnetic waves traveling in the z-direction, with the electric field vector pointing in 
the ^/-direction, thus representing a two-dimensional wave propagation problem. 

The forward problem, consisting of obtaining reflectance data for a known structure, may 
be formulated in terms of transfer matrices, providing reflection and transmission matrices 
for each layer as well as the whole structure. The corresponding inverse problem, consisting 
of recovering the properties of the unknown structure from the reflectance data, requires 
data for several wave numbers and a range of frequencies. 

The key step of the synthetic layer-stripping algorithm is to combine reflectance data for 
different frequencies so as to synthesize data from a short pulse at the time the pulse hits 
the surface. The data may then be applied to reconstruct the refractive index profile in the 
leftmost layer by identifying the local, time-domain Fresnel reflection at each point. Once 
the first layer is known, its transfer matrix may be computed, and the impact of the layer 
on the reflectance data eliminated. The same idea is then applied to the modified reflectance 
data, and continuing in the same way we are, at least in principle, able to reconstruct the 
whole structure. 
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II. The forward problem 

We are considering photonic crystals consisting of homogeneous, dielectric and nonmagnetic 
(/i = /io) materials, where the dielectric permittivity, e{x,z), is varying in the x- and z- 
directions, and being constant in the y-direction. The structure consists of layers orthogonal 
to the 2;-axis, with e constant with respect to z (and y) in each layer. We assume that the 
structure contains layers of finite thickness occupying the space from zo = to z^. Layer 
i spans < z < Zi, with a thickness Aj = — and permittivity ej(x). Outside the 
structure we assume that it is vacuum. By inserting the Fourier transform in time of the 
electric and magnetic fields, e.g. 



1 r 



e(x,z,uj)e-"^'duj, (1] 



and similarly for H, into Maxwell's equations, we obtain 

V X e{x, z,u) — iufioh{x, z,u) = 0, (2) 

V X h{x, z,u) + iujeoe{x, z)e{x, z,uj) = 0, (3) 

V ■ {e{x,z)e{x,z,u)) = 0, (4) 

V-h{x,z,uj) = 0. (5) 

In addition, from the standard continuity conditions for E and H, parallel (||) and orthogonal 
(±) to the layer boundaries, the following conditions have to apply 



eiei 


= e2e: 


ej 


- e'l 


hi 


= h2. 
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(6) 
(7) 



Omitting the dependence of u, direct substitution shows that the Fourier transformed 
Maxwell's equations in the present case admit solutions of the form 

e{x,z) = e{x,z)y, (9) 
w N 1 f de{x,z)^ de(x,z)^\ 

h{x,z) = -. ^^x+ ^ ' , 10 

zujfiQ \ oz ox J 

provided e (x, z) satisfies the scalar Helmholtz equation, 

+ |^ + 6(x,z)A;2e = O, (11) 

where k = u / Cq, = l//ioeo, and e as as well as its partial derivatives are continuous across 
the layer boundaries. 
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The structure will be probed using the plane polarized waves in Eqs. ([9]) and (fTOj) . and all 
transmitted and reflected waves will have the same polarization. 

The direct problem consists of solving Eq. ffTTj) with appropriate boundary conditions. 
These will be plane waves entering from the half space z < (thus, having positive z- 
component wave numbers), whereas no waves are supposed to enter the structure from the 
region z > Zn- 

III. Solution by transfer matrices 

Let L be the extension of the computational domain in the x-direction. By imposing L- 
periodic boundary conditions, we may expand the electric field and the permittivity into the 
Fourier series 



e{x,z) = e{x,z)y 



m). 



e{x,z) = ^e(™n^)exp(i4'")a;). 



(12) 
(13) 



where /ci"^^ = 2'Km/L. The Fourier components are given as 



e^'^Hz) 



e[x, z) exp(— z 



k^^^x)d 



X. 



By inserting the Fourier series into Helmholtz' equation we obtain 

(A:i'"))'E(™)(2) + A;2^e("')(2)E(™')(^) = 0, 



J2 77'{m)/' \ 



dz^ 



(14) 



(15) 



for each m. Let E = {E^'^^}jn£Z and use = ^/ci™''^ + ^A;!'"^^ . We may then write Eq. 
f|T5|) as the matrix equation 



d'E(^) 
dz^ 



{kl + Viz)) E{z) = 0, 



(16) 



where = diag(/ci'"'*), and V(z) is the infinite dimensional Toeplitz matrix operator defined 



as 



v(^) = -kH + e 



e{o) 



e{-2) 

e{0) 



.(2) .(1) .(0) 



(17) 
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Eq. ( IT6|) may be decomposed into a first order system by writing E{z) = £+(2) + E (z) 
and require 

dE+r^) 



dz 
dE-{z) 
dz 



tKB+iz) + 2(2k,)-iV(2)(E+(^) + E-(^)), 
-zKE^iz) - z{2k,y'V{z){E+{z) + E"(z)). 



18a 



(18b) 



Eq. (fT6l) follows easily from Eqs. f ll8ap and ( 118bl) after summation and differentiation. Note 
that outside the structure V equals zero, and the solution to Eq. f|T8l) are forward and 
backward going waves. 
By writing 

^ tK + ti2K)-'V{z) i{2K)-'V{z) 
-i{2k,y^V{z) -ik, - i{2k,y^\{z) 

Eq. f[T5|) can be brought into the compact matrix form 

d^(z) 









, C{z) = 







(19) 



dz 



C(z)^(z). 



(20) 



Since the permittivity is independent of z within a layer, the matrix C will be constant for 
each layer. Thus, Eq. f l2Q|) can be integrated to 

^{zb) = exp[{zb-Za)C]^iza) (21) 

for Za and zi, inside the same layer. Let Cj be the C-matrix for layer i. By defining 

M, = exp(A,Ci), (22) 

we obtain 

*(z,) = M,*(z,_i). (23) 

Since e and de/dz are continuous across the layer boundaries, E+ and E~ also must be 
continuous here. By successive apphcations of Eq. (123!) we are able to propagate through the 
structure from zq to z^, 



*(^^) = J]M,*(^o) = M*(^o). 

i=N 

The boundary conditions at z = zq = are contained in 



E+(zo) 
E-(zo) 



(24) 



(25) 



where E~^{zo) is given by the probing waves, and E (zq) is unknown. At the other end of 
the structure, 

' E^{zn) 




(26) 
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where E'^{zn) has to be determined. If we partition the matrix M according to the definition 
of*. 

Mil Mi2 
M21 M22 



M 



(27) 



and consider the matrix equation 



" E+(^^) " 




" Mil 


M12 




" E+(^o) " 







_ M21 


M22 




_ E-(^o) _ 



(28) 



the formal solution is easily seen to be 



E-(zo) = RE+(zo), (29) 
E+(^;v) = TE+(^o), (30) 

where R and T are the so-called reflection and transmission matrices for the structure, 

R = -M^2'M2i, (31) 
T = Mil - Mi2M22^M2i. (32) 

The refiection and transmission matrices provide 'E~(zq) and E'^{zi^) for all possible input 
E+(zo). 

To solve the inverse problem, we need refiectance data for all possible incident wave num- 
bers. Such a set of excitation-response pairs can be described by the following equation 



T 




I 




= M 









R 



(33) 



Here, each column z of I corresponds to an experiment where the incident field amplitude 
is 1 for one of the Fourier components and zero for the others. The ith column of R is the 
refiection at zq, and the ith column of T is the corresponding transmission at zj^. 

IV. The inverse problem 

The inverse problem is solved by combining refiectance data for a range of frequencies, 
[^1,1^2], so as to synthesize data from a short pulse at the time the pulse hits the surface. 
The permittivity of the layer may then be recovered as described below if the pulse width in 
time, O (10 {u2 — oui) ), is shorter than the round-trip travel time in the layer, 2Aj/c. In 
this study we shall, for simplicity, assume that the layer thicknesses, Aj, are known. It is, in 
principle, possible to do the layer-stripping without knowledge of the layer thickness. This 
will be briefiy discussed this in Sec. IVII From the permittivity and the thickness it is now 
possible to compute the transfer matrix for the layer, and the impact of the layer removed 
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from the reflection data. Thus, we obtain reflectance data for the same structure, but without 
the leftmost layer, and may then repeat the steps above until the entire structure has been 
reconstructed. 



IV. A. Computing the permittivity 

Let us assume that we have reflectance data from an experiment where the incident wave is 
a plane wave pulse, F(t — z/cq), see Fig. [1] The incident wave enters from z < 0, and first 
hits the leftmost layer of the structure. The wave speed in the first layer is Ci{x). 

Let us consider a small neighborhood around the point x, and assume that the permittivity 
in the leftmost layer is varying slowly enough, so that we may take it to be constant and 
equal to ei(x) inside the area we consider. Furthermore, we restrict our calculation to a small 
time interval so that the reflected and transmitted waves are only affected by the first layer. 

At the left side of the boundary, the electric field will be a sum of the incoming and 
the reflected wave, and at the right side, the field consists of the transmitted wave. Thus, 
continuity of the transversal electromagnetic fields gives us 

F{t --) + R,{x)F{t +-) = T,{x)F{t - ^), (34) 

Co Co Ci{X) 

and 

--F\t - ^) + ^F'it + ^) = -^^F'it - 4-), (35) 
Co Co Co Co ci{x) ci{x) 

for some reflection coefficient Ri{x). Obviously, the incident wave pulse F needs to be short 
for Eqs. (!34|) and (135|) to hold. For z = and < t <^ Ai/cq we thus have from Eqs. (l34l) 
and (J35]) 

1 + Ri{x) = Ti(x), (36) 
-1 + Riix) = -Ti(x)^, (37) 

Ci(x) 

leading to the following expression for the permittivity in the leftmost layer: 

Eq. (l38l) can be interpreted as a local Fresnel equation; it connects the local reflection coef- 
ficient Ri{x) and the local permittivity ei{x). 

IV. B. Layer-stripping 

For the layer-stripping we shall assume that the layer thickness is known. Once we have 
computed ei(x), we obtain the transfer matrix Mi from Eq. (l22l) . As in Eq. fl23|) . the forward 
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and backward going waves just before the second layer can be computed as 



" E+(zO " 


= Ml 






= Ml 


I 




E- 




R 



Note that here E+ and E are matrices as in Eq. 
structure, Eqs. (129|) and (130|) . we may define 



(39) 



. Similar to the solution for the full 



E-(^i) 



TE+(^i) 
RE+(^i^ 



and thus the new refiection matrix can be found as 

R = E-(^i)E+(^i)-i. 



(40) 
(41) 



(42) 



V. Algorithm 

For the numerical experiments, which we will discuss in Sec. \V1\ we first solve the forward 
problem to obtain valid refiectance data Ti{k^\k^\uj) for each of A*"^ frequencies in the 
span [ci;i,ci;2]. Here /ci*" denote the incident wave numbers, and k^'^ the refiected wave 
numbers. 

For the inverse problem, the first step is to synthesize a time localized pulse in space and 
time to be applied in Eq. ( I38l) . Only data for incident waves normal to the structure are used 
in the identification of the layers. We start with a transformation from ki"^^ to x. 



r{x,u) = Y,mt"> = 0,k^r\^)e 



(43) 



where 



Ay™ 



27rm 



m E Z, and u G [ci;i,Ci;2]. 



(44) 



Ideally, one would now carry out an inverse Fourier transform in order to obtain the refiected 
field in the space- and time-domains resulting from an incident delta-pulse. 



R{x,t) 



1 

2^ 



r X, u e 



-itut 



du. 



(45) 



However, since we only know r{x,u) in the frequency interval [co'i,Co'2]; we need to choose 
a window function W{u) with support in that interval. The corresponding pulse in time is 
denoted w{t), and the experience with some standard window functions are discussed in the 
next section. The synthetic response will then be 



1 

R^{x,t) = — r{x,u)W{u)e-''^'du. 
27r J^, 



(46) 
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For z = and < t ^ Ai, we see from Eq. that {x,t) = Ri (x) w (t). Hence, for 
t = 0, we may write 

_ R^ix,0) _ j::rix,u;)W{u)du 

w{0) ~ JW{u)du • ^ 

The permittivity e (x) is then obtained from Eq. ( l38l) . The layer-stripping has been discussed 
in Sec. irvTBl 

The numerical method may now be summarized in a pseudo code, as follows. 

Input to the algorithm: Reflection matrices R(/ci"^ ■* = 0, k^\u), for u G [aJi,a;2]. 

1 Transform reflection matrices to the space/frequency domain using Eq. fH51) . Next, 
determine the reflection coefficient at t = using Eqs. (H6|) - (147|) . The permittivity 
now follows from Eq. (138|) . 

2 Use e(x) and Eq. (fT4|) to find the Fourier components. Then the transfer matrix can 
be computed from Eqs. (fTTD. ffTOj) and (l22|) . 

3 Knowing the transfer matrix, the forward and backward going waves before the next 
layer can be computed from Eq. (1391) . 

4 Compute the new reflection matrix, using Eq. (142|) . 

5 Return to step 1, to compute the permittivity in the next layer, until all layers have 
been found. 

VI. Numerical experiments 

For simplicity, we scale eo and /iq (in free space) to 1, and hence cq = 1. Now, the dispersion 
relation in free space reduces to k = u, and the wavelengths will be A = 2ti/uj. The length 
of the synthesized pulse will depend on the frequency interval [ci;i,ci;2], and be of order 
27r {uj2 — ooi)~^ ■ One must require that the pulse length is shorter than the round-trip travel 
time in the first layer, i.e. 

< (48) 



U2-UJ1 ~ Cl 

Thus, to be able to reconstruct the structure, we need to make sure that the frequency band 
is wide enough. Through simulations we have seen that for smaller frequency bands, the 
error fluctuates as the bandwidth grows. However, as soon as the frequency span is wide 
enough, the calculations are stable. 

As mentioned in Sec. |Vl we use Eq. fl33|) to compute reflection matrices. To compute 
the reflection matrix and solve the inverse problem, the number of Fourier components, 
M, the computational domain, L, and the A'^ frequencies in the interval [a;i,a;2], must be 
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chosen. Now, the reflection and transmission matrices will have dimension M x M and the 
propagation matrix M dimension 2M x 2M. Note that given the resolution in the x-direction 

Ax = ^, (49) 
the maximum Nyquist wave number is given as 

maxl, = ^. (50) 
L 

Since 

K = ^fc^^ = ^/^^^^x. (51) 

evanescent modes in vacuum occur when maxA;^; > u. 

As we see in Eq. ( H2l) . the updated reflection matrix is a product of the backward traveling 
waves and the inverse of the forward traveling waves. In a homogeneous medium, the forward 
and backward traveling waves have the z-dependence exp(z/c2z) and exp(— z/c^z), respectively. 
Thus, when we have evanescent modes with kz = iK for some > 0, we expect a noise 
amplification factor of the order of exp(2ii'A) when one layer is stripped off. Because of this, 
we restrict our frequency band to avoid a too large K. 

Note that if the variations in the x-direction are fast, we will need a finer resolution in x. 
To obtain this, we must either make M larger, or L smaller, which again implies that the 
maximum wave number, max/c^;, increases. The result is that we need higher frequencies to 
avoid evanescent modes. Thus, it is convenient to restrict the attention to structures that 
vary slowly with respect to x. We will, however, in Example 2 see that we are also able to 
reconstruct a structure with fast variations. 

There are three effects that mainly contribute to the limitations in this method. The first 
one is the effect of evanescent modes, which was discussed above. The second is related to the 
contrast between the minimum and maximum refractive indices, rii and n2. As the contrast 
grow higher, more of the light is reflected by the flrst layers. This will in turn make the 
calculations less accurate [8l|9] . The last effect is also related to the bandwidth. As discussed 
in the beginning of this section, we need the pulse length to be smaller than the layer round- 
trip travel time, thus imposing a lower limit to the bandwidth. Conversely, if the bandwidth 
is flxed, the layer thicknesses have to be large enough. 

As discussed in Sec. |Vl it is necessary to shape the reflection data using some window 
function. Numerical experiments applying three different window functions on the interval 
[cji,cj2] are presented in Sec. IVI.AI As expected, the experiments favor smoother windows 
over the simple rectangular window. Consequently, a Hanning window has been applied for 
the rest of the computations. 

In all simulations we have assumed that we know the layer thickness. It would, at least 
in principle, be possible to reconstruct the layer thickness. To do this, one need to use a 
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fictitious layer thickness, A2;, which must be small compared to the expected layer thickness. 
Then this Az can be used in the computations, to reconstruct each layer piece by piece. 

All the calculations where done using MATLAB on an Intel Core 2 Quad 2.83 GHz com- 
puter. For a typical calculation, the runtime was approximately 260 s. About 180 s was used 
on calculating the forward problem, i.e. the reflection matrices R(A;i™ \ k^\u), the calcula- 
tion of the permittivity is neglectable, and the calculation and removal of the leftmost layer 
took about 80 s. 

VI. A. Example 1 

In this example we have chosen a structure where the permittivity within the layers is given 

as 

X 

e(x) =r/±7cos(-), (52) 

and where the two different layers are obtained by ahernating between + and — . The layer 
thickness has been set to 7r/2, and the parameters r] and 7 are chosen so that 

mine = n\ = 1.0, (53) 
maxe = (54) 

i.e. rii is kept constant, while n\ is varied. For all the experiments presented here, M = 300 
and L = 100, giving a resolution in a; of 1/3 and a maximum kx of Sir. The choice of frequency 
band depends on whether we want to include evanescent modes in the computations, how 
large K we can tolerate, and the layer thickness. 

We first consider the choice of frequency window functions, Eq. f H6|) . by testing three dif- 
ferent functions, a simple rectangular window, the Banning window and the Tukey (tapered 
cosine) window. The comparisons were carried out for two different values oi n\ = 1.2 
for the first test, and n\ = 2.0 for the second. The real difference in the choice of window 
functions can only be seen in the reconstruction of the second layer. Therefore, the results, 
displayed in Fig. [21 are only shown for this layer. The rectangular window give poor results 
already for the lowest refractive index, whereas the two other choices give a nice recon- 
struction. For the higher index, we see that also the Tucky window is beginning to give less 
accurate results, while for the Banning window there is almost no visible difference from the 
exact curve. For the rest of the numerical examples, we have therefore chosen to use the 
Banning window. 

Let us now turn to the effect of including evanescent modes. In these experiments, the 
permittivity was first chosen such that = 1.0 and n"^ = 1.05, and secondly we had nl = 1.0 
and rig = 2.0. The frequency bandwidth was set constant to 10, but the lower and upper 
frequencies were changed to include none, or some evanescent modes. Again, the first layer 
is nicely reconstructed, so we only show results for the second layer in Fig. |3l Including some 
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evanescent modes improves the results, but if K gets too large, there are fluctuations in the 
reconstruction. Note that as n2 gets higher, a larger K is acceptable . This can be explained 
by the fact that the fleld now is locally non-evanescent, since a local in the medium would 



Bigger contrasts in the refractive index may be simulated by keeping the lower index 
constant at n\ = 1.0, while the upper index is changed. The results are shown in Fig. HJand 
we see that the flrst layer is well reconstructed for all choices of n"^, while some fluctuations 
in the second layer are visible when n| > 2. 

Finally, we demonstrate that it is possible to reconstruct more than two layers. Fig. [5] 
shows the error, given as |ecomp — e| for each layer in a structure consisting of four layers. 
The permittivity was chosen between nl = 1.0 and = 1.05. 

VLB. Example 2, a square function 

In this example we have chosen the permittivity to be a square function with period 27r in 
each layer. To get two different layers, and a permittivity between rii and 77-2, we let 



This example has been chosen to show that we are able to reconstruct a less smooth func- 
tion. The fact that the permittivity changes fast seems to be in conflict with the discussion 
in Sec. IVIt but as we see from Fig. [6] we actually manage to reconstruct two layers of this 
structures as well. The reason it work is that we only use the low frequency components of 
the square. The results can be seen in Fig. El 

VII. Conclusion 

In this paper we have presented a method for reconstructing 2d photonic structures layer- 
by-layer. In principle, the reconstruction is exact, and we have shown through numerical 
examples that we are able to reconstruct different structures consisting of several layers. 

There are essentially three, fundamental mechanisms that limit the accuracy in practice: 
The presence of evanescent modes, accumulated reflection, and limited bandwidth. The fact 
that evanescent modes lead to inaccuracies means that either the probing frequencies must be 
sufficiently high, or the spatial transversal frequencies of the structure must be sufficiently 
small. The second limitation is a result of the fact that if the transmission through the 
structure is too small, little light reaches the back end. Then the back end has little influence 



be h = {n^co^ - klY'\ 




(55) 



(56) 
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on the reflection data, and cannot be reconstructed accurately. Finally, the bandwidth must 
be sufficiently large such that the synthetic, incident pulse is shorter than the round-trip 
time in all layers. Alternatively, for a fixed, available bandwidth, the structure must vary 
sufficiently slowly in the longitudinal direction z. In the numerical examples, for a center 
wavelength 1 |j,m the normalized frequency interval [9, 19] corresponds to the wavelength 
interval [0.74, 1.56] |J.m. The layer thicknesses in the example become 3.5 |a.m. If the available 
bandwidth is reduced to 100 nm, the layers must be 5-10 times thicker to achieve the same 
accuracy. 

References 

1. E. Yablonovitch, "Inhibited Spontaneous Emission in Solid-State Physics and Electron- 
ics", Phys. Rev. Lett. 58, 2059-2062 (1987). 

2. S. John, "Strong Localization of Photons in Certain Disordered Dielectric Superlattices" , 
Phys. Rev. Lett. 58, 2486-2489 (1987). 

3. J. D. Joannopoulos, R. D. Meade and J. N. Winn, Photonic Crystals: Molding the Flow 
of Light (Princeton University Press, 1995) 

4. A. M. Bruckstein, B. C. Levy and T. Kailath, "Differential Methods in Inverse Scatte- 
ring", SIAM J. Appl. Math. 45, 312-335 (1985). 

5. A. E. Yagle and B. C. Levy, "Layer-stripping solutions of multidimensional scattering 
problems", J. Math. Phys. 27, 1701-1710 (1986). 

6. J. Skaar, L. Wang and T. Erdogan, "On the Synthesis of Fiber Bragg Gratings by Layer 
Peeling", IEEE J. Quantum Electron. 37, 165-173 (2001). 

7. O. H. Waagaard and J. Skaar, "Inverse Scattering in Multimode Structures", SIAM J. 
Appl. Math. 68, 311-333 (2007). 

8. A. M. Bruckstein, I. Koltracht, and T. Kailath, "Inverse scattering with noisy data", 
SIAM J. Sci. Stat. Comput. 7, 1331-1349 (1986). 

9. J. Skaar and R. Feced, "Reconstruction of gratings from noisy reffection data", J. Opt. 
Soc. Am. A 19, 2229-2237 (2002). 



13 




a) t < h)t = 0+ 

Fig. 1. (a) A plane wave pulse is incident to a layered 2d structure. For the 
structure different colors indicate different refractive indices, (b) Immediately 
after t = the field has only been affected by the first layer; thus we may 
identify the first layer from the first part of the reflected field in the time- 
domain. AKSFigl.eps. 
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Fig. 2. Reconstruction of the second layer for different window functions, in (a) 
nl = 1.2 and in (b) = 2.0. Other parameters: M = 300, L = 100, = 100, 
u G [9,19], Az = 7r/2, nf = 1.0. Note that only parts of the computational 
domain L are shown in the plots. AKSFig2.eps. 
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Fig. 3. Reconstruction of the second layer for different ui and 002, but with 
constant bandwidth. In (a) n\ = 1.05, and in (b) nl = 2.0. Other parameters: 
M = 300, L = 100, = 100, Az = tt/2, n\ = 1.0. AKSFig3.eps. 
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Fig. 4. Reconstruction of the first layer, (a), and the second layer, (b), for 
different choices of n2- Other parameters: M = 300, L = 100, u G [9,19], 
= 100, Az = 7r/2, nf = 1.0. The solid black lines represent the exact 
permittivity. AKSFig4.eps. 
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Fig. 5. The error in each layer, given as |ecomp — e|. Parameters: M = 300, 
L = 100, cu e [9, 19], A^^ = 100, Az = 7r/2, nj = 1.0, nj = 1.05. AKSFig5.eps. 
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Fig. 6. Reconstruction of the square function for the first layer, (a), and for 
the second layer, (b). The computations were done for different ui and U2, but 
with constant bandwidth. Other parameters: M = 300, L = 100, N^^ = 100, 
Az = 7r/2, nj = 1.0, nl = 1.05. AKSFig6.eps. 
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